<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Adding more complex constraints on a model </title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S4 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S14 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S16 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Adding more complex constraints on a model </span></h1><h2  class = 'S1'><span>Authors: Thomas Pfau,Life Sciences Reasearch Unit, University of Luxembourg, Luxembourg.</span></h2><h2  class = 'S1'><span>Reviewers: </span></h2><h2  class = 'S2'><span>INTRODUCTION</span></h2><div  class = 'S3'><span>The COBRA Toolbox offers the possibility to add additional constraints to a model that are not direct flux constraints. One example is already indicated in the constraining models tutorial which introduces coupling constraints[1], but there are plenty of other possibilities that can be achieved using constraints.</span></div><div  class = 'S3'><span>The tools introduced in this tutorial are mainly aimed at developers who want to implement complex algorithms or formalism within the COBRA Toolbox environment, but the basics are also useful to integrate specific properties of a model or specific literature data. The section "The COBRA Model structure" is aimed at developers and can be skipped by users who don't to go into the details.</span></div><div  class = 'S3'><span></span></div><h2  class = 'S2'><span>The COBRA Model structure</span></h2><div  class = 'S3'><span>In general, the COBRA toolbox represents a model in a struct using multiple fields to represent the different properties of the model. Of particular interest for analysis is the stoichiometric matrix S (which represents the metabolic stoichiometries of the reactions in the model). Under normal circumstances, this matrix can directly be translated into the linear problem matrix A of a problem struct as used by the Toolbox. However, since there are several algorithms which aim at manipulating reactions or metabolites, it is important to draw a distinction between the stoichiometric matrix S and any other constraints or variables which can be added to the model.</span></div><div  class = 'S3'><span>Therefore 3 additional matrices exist, where such information can be stored:</span></div><div  class = 'S3'><span>E - The matrix indicating the influence of additional variables on metabolite levels.</span></div><div  class = 'S3'><span>C - The matrix of additional constraints on the reactions, which are not derived from the steady state condition and</span></div><div  class = 'S3'><span>D - The matrix of interactions between additional Constraints and additional variables.</span></div><div  class = 'S3'><span>Overall A COBRA model will be translated into a Linear problem by combining these matrices in the following way:</span></div><div  class = 'S4'><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;block&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;LPproblem&lt;/mi&gt;&lt;mo&gt;.&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;A&lt;/mi&gt;&lt;mo&gt;=&lt;/mo&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;mrow&gt;&lt;mo&gt;[&lt;/mo&gt;&lt;mtable&gt;&lt;mtr&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;model&lt;/mi&gt;&lt;mo&gt;.&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;S&lt;/mi&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;model&lt;/mi&gt;&lt;mo&gt;.&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;E&lt;/mi&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;/mtr&gt;&lt;mtr&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;model&lt;/mi&gt;&lt;mo&gt;.&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;C&lt;/mi&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;mtd&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;model&lt;/mi&gt;&lt;mo&gt;.&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;D&lt;/mi&gt;&lt;/mrow&gt;&lt;/mtd&gt;&lt;/mtr&gt;&lt;/mtable&gt;&lt;mo&gt;]&lt;/mo&gt;&lt;/mrow&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-14px"><img src="" width="247.5" height="39" /></span><span> </span></div><div  class = 'S3'><span>The Cobra toolbox allows one sided inequality constraints and equality constraints on a model. i.e. there is currently no mechanism to add a constraint like:</span></div><div  class = 'S4'><span texencoding="c_{low} \le v_{R1} \le c_{high}" style="vertical-align:-6px"><img src="" width="104" height="20" /></span></div><div  class = 'S3'><span>directly. Instead it would be necessary to add the two constraints:</span></div><div  class = 'S4'><span texencoding="c_{low} \le v_{R1}" style="vertical-align:-6px"><img src="" width="61" height="20" /></span></div><div  class = 'S3'><span>and</span></div><div  class = 'S4'><span texencoding=" v_{R1} \le c_{high}" style="vertical-align:-6px"><img src="" width="64.5" height="20" /></span></div><div  class = 'S3'><span>which is commonly what solvers will translate these constraints into anyways. </span></div><div  class = 'S3'><span>It is best practice to not manually alter the sizes of these matrices, but to use the appropriate model manipulation functions (addReaction/addMultipleReactions, addMetabolite/addMultipleMetabolites, addCOBRAVariables, addCOBRAConstraints). These functions will ensure, that the model structure stays in sync, and all necessary fields are updated appropriately. There is a convenience function to create the corresponding LPproblem struct to the given model: buildOptProblemFromModel(). This function will create and initialize all necessary fields and builds the corresponding LP problem. Any additional modifications of the LP problem imposed by an algorithm should be done on this LP and not on the original model structure.</span></div><div  class = 'S3'><span></span></div><div  class = 'S3'><span></span></div><h2  class = 'S2'><span>PROCEDURE</span></h2><div  class = 'S3'><span>Initially, we will load a model that we want to modify by adding a few additional constraints. The model used will be the simple E.Coli core model:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: normal"><span >initCobraToolbox(false)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="0BE58C15" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="829" data-hashorizontaloverflow="true" data-testid="output_0" style="max-height: 261px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">      _____   _____   _____   _____     _____     |
     /  ___| /  _  \ |  _  \ |  _  \   / ___ \    |   COnstraint-Based Reconstruction and Analysis
     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - 2018
     | |     | | | | |  _  { |  _  /  |  ___  |   |
     | |___  | |_| | | |_| | | | \ \  | |   | |   |   Documentation:
     \_____| \_____/ |_____/ |_|  \_\ |_|   |_|   |   <a href="http://opencobra.github.io/cobratoolbox" style="white-space: normal; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">http://opencobra.github.io/cobratoolbox</a>
                                                  | 

 &gt; Checking if git is installed ...  Done (version: 2.7.4).
 &gt; Checking if the repository is tracked using git ...  Done.
 &gt; Checking if curl is installed ...  Done.
 &gt; Checking if remote can be reached ...  Done.
 &gt; Initializing and updating submodules (this may take a while)... Done.
 &gt; Adding all the files of The COBRA Toolbox ...  Done.
 &gt; Define CB map output... set to svg.
 &gt; TranslateSBML is installed and working properly.
 &gt; Configuring solver environment variables ...
   - [*---] ILOG_CPLEX_PATH: /opt/ibm/ILOG/CPLEX_Studio128/cplex/matlab/x86-64_linux
   - [*---] GUROBI_PATH: /opt/gurobi800/linux64/matlab
   - [----] TOMLAB_PATH: --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: normal; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   - [-*--] MOSEK_PATH: /home/thomas/mosek/8
   Done.
 &gt; Checking available solvers and solver interfaces ... Done.
 &gt; Setting default solvers ... Done.
 &gt; Saving the MATLAB path ... Done.
   - The MATLAB path was saved as ~/pathdef.m.

 &gt; Summary of available solvers and solver interfaces

			Support 	   LP 	 MILP 	   QP 	 MIQP 	  NLP
	----------------------------------------------------------------------
	gurobi       	active        	    1 	    1 	    1 	    1 	    -
	ibm_cplex    	active        	    1 	    1 	    1 	    - 	    -
	tomlab_cplex 	active        	    0 	    0 	    0 	    0 	    -
	glpk         	active        	    1 	    1 	    - 	    - 	    -
	mosek        	active        	    1 	    - 	    1 	    - 	    -
	matlab       	active        	    1 	    - 	    - 	    - 	    1
	cplex_direct 	active        	    0 	    0 	    0 	    0 	    -
	dqqMinos     	active        	    1 	    - 	    - 	    - 	    -
	pdco         	active        	    1 	    - 	    1 	    - 	    -
	quadMinos    	active        	    1 	    - 	    - 	    - 	    -
	qpng         	passive       	    - 	    - 	    1 	    - 	    -
	tomlab_snopt 	passive       	    - 	    - 	    - 	    - 	    0
	lp_solve     	legacy        	    1 	    - 	    - 	    - 	    -
	----------------------------------------------------------------------
	Total        	-             	    9 	    3 	    5 	    1 	    1

 + Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.


 &gt; You can solve LP problems using: 'gurobi' - 'ibm_cplex' - 'glpk' - 'mosek' - 'matlab' - 'dqqMinos' - 'pdco' - 'quadMinos' 
 &gt; You can solve MILP problems using: 'gurobi' - 'ibm_cplex' - 'glpk' 
 &gt; You can solve QP problems using: 'gurobi' - 'ibm_cplex' - 'mosek' - 'pdco' - 'qpng' 
 &gt; You can solve MIQP problems using: 'gurobi' 
 &gt; You can solve NLP problems using: 'matlab' 

 &gt; Checking for available updates ...
 --&gt; You cannot update your fork using updateCobraToolbox(). [b08fe8 @ ConstraintFixes].
     Please use the MATLAB.devTools (<a href="https://github.com/opencobra/MATLAB.devTools" style="white-space: normal; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">https://github.com/opencobra/MATLAB.devTools</a>) to update your fork.</div></div></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >model = getDistributedModel(</span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%Create a copy for comparisons</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model_orig = model;</span></span></div></div></div><div  class = 'S10'><span>We will add a restriction on the activity of the two aconitase proteins present in E.Coli (aconA and aconB), which catalyse two steps in the citric acid cycle (see below)</span><img class = "imageNode" src = "" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S3'><span></span></div><div  class = 'S3'><span>From UNIPROT[2], the Km value of aconitase acting on citrate is 6.13 and 23.8 umol/min/mg for aconA and aconB respectively. For cis-aconitate the values are 14.5 and 39.1 umol/min/mg, respectively. Since most models assume fluxes having the unit mmol/gDW/hr we will have to convert the units to mmol/hr instead of umol/min:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >aconACit = 6.13 / 1000 * 60;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconAAcon = 14.5 / 1000 * 60;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconBCit = 23.8 / 1000 * 60;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >aconBAcon = 39.1 / 1000 * 60;</span></span></div></div></div><div  class = 'S10'><span>From Wiśniewski and Rakus[3] the amount of aconA and aconB per mg E.Coli sample is ~4.05 pmol/mg and 95.95 pmol/mg respectively, with a weight of 97.676 kDa and 93.497 kDa, respectively.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >aconAmol_per_g = 4.05 * 1000 * 1e-12;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconBmol_per_g = 95.95 * 1000 * 1e-12;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconA_molWeight = 97.676 * 1e3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconB_molWeight = 93.497 * 1e3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >aconAAmount = aconAmol_per_g * aconA_molWeight / 0.3; </span><span style="color: rgb(2, 128, 9);">% divided by 0.3 to account for </span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >                                                      </span><span style="color: rgb(2, 128, 9);">% the non water fraction assuming</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >                                                      </span><span style="color: rgb(2, 128, 9);">% 70% water.</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >aconBAmount = aconBmol_per_g * aconB_molWeight / 0.3;</span></span></div></div></div><div  class = 'S10'><span>Now, there are two genes which code for aconitase in the ecoli core model: b0118 (aconB) and b1276 (aconA). Both the citrate hydratase and the aconitatde dehydratase have the same GPR rule: (b0118 or b1276), so they can both use either enzyme. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >aconAgene = </span><span style="color: rgb(170, 4, 249);">'b1276'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >aconBgene = </span><span style="color: rgb(170, 4, 249);">'b0118'</span><span >;</span></span></div></div></div><div  class = 'S10'><span>We would like to add a constraint that not only restricts the activity of these two reactions, but also ensures that the turnover rates are considered. One mg of aconA can support a total flux of ~0.367 mmol/h of ACONTa</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: normal"><span >printRxnFormula(model,</span><span style="color: rgb(170, 4, 249);">'rxnAbbrList'</span><span >,{</span><span style="color: rgb(170, 4, 249);">'ACONTa'</span><span >},</span><span style="color: rgb(170, 4, 249);">'gprFlag'</span><span >, true);</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="A041C606" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="18" data-hashorizontaloverflow="false" data-testid="output_1" style="max-height: 261px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ACONTa	cit[c] 	&lt;=&gt;	acon-C[c] + h2o[c] 	(b0118 or b1276)</div></div></div></div></div><div  class = 'S10'><span>OR .087 mmol/hr through ACONTb.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: normal"><span >printRxnFormula(model,</span><span style="color: rgb(170, 4, 249);">'rxnAbbrList'</span><span >,{</span><span style="color: rgb(170, 4, 249);">'ACONTb'</span><span >},</span><span style="color: rgb(170, 4, 249);">'gprFlag'</span><span >, true);</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="728D4AA0" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="18" data-hashorizontaloverflow="true" data-testid="output_2" style="max-height: 261px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ACONTb	acon-C[c] + h2o[c] 	&lt;=&gt;	icit[c] 	(b0118 or b1276)</div></div></div></div></div><div  class = 'S10'><span>However, it will not be able to do both at the same time. Therefore, we can assume that in addition to its normal metabolites, the reaction also consumes some of the enzyme (for this time step). However, both reactions can also be catalysed by aconB, so they might actually not use any of aconA. Therefore, we need additional constraints that represent the activity through these two reactions. In addition, we need variables that represent the efficiency of one mg of protein to catalyse the reactions.</span></div><div  class = 'S3'><span>Essentially, we have to encode:</span></div><div  class = 'S3'><span>A: The availability of the given proteins </span></div><div  class = 'S3'><span>B: The efficiency of a protein to catalyse a specific reaction</span></div><div  class = 'S3'><span>C: The usage of that protein by the respective reaction.</span></div><div  class = 'S3'><span>And we will do so in this tutorial</span></div><h2  class = 'S2'><span>Adding availability variables.</span></h2><div  class = 'S3'><span>The rxns field is intended to only represent reactions from the model, so if an additional variable is required for a specific task, this variable should not be generated in the rxns field, but in a distinct field for this kind of variables. Therefore, we will add availability variables, which can be thought of as exchange reaction for the enzymes. To do so, we use the following function</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%addCOBRAVariables(model, idList, varargin)</span></span></div></div></div><div  class = 'S10'><span>The function takes an existing model, and a set of variable ids (</span><span style=' font-weight: bold;'>idList</span><span>) and generates all required fields for these variables. You can modify the properties of the variables by the following parameters which can either be supplied as parameter/value pairs:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >addCOBRAVariables(model,{</span><span style="color: rgb(170, 4, 249);">'newVar'</span><span >},</span><span style="color: rgb(170, 4, 249);">'lb'</span><span >, -5,</span><span style="color: rgb(170, 4, 249);">'ub'</span><span >,3);</span></span></div></div></div><div  class = 'S10'><span> or as a parameter struct:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >params = struct();</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >params.lb = -5;</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >params.ub = 3;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >addCOBRAVariables(model,{</span><span style="color: rgb(170, 4, 249);">'newVar'</span><span >},params);</span></span></div></div></div><div  class = 'S10'><span>The available parameters are:</span></div><ul  class = 'S13'><li  class = 'S14'><span style=' font-weight: bold;'>'lb'</span><span> - The lower bound(s) of the variable(s) (default: -1000)</span></li><li  class = 'S14'><span style=' font-weight: bold;'>'ub'</span><span> - The upper bound(s) of the variable(s) (default: 1000)</span></li><li  class = 'S14'><span style=' font-weight: bold;'>'c'</span><span> - The objective coefficient(s) of the variable(s) (default: 0)</span></li><li  class = 'S14'><span style=' font-weight: bold;'>'Names'</span><span> - Descriptive name(s) of the variable(s) (default: idList)</span></li></ul><div  class = 'S3'><span>We will now add the variables 'aconA' and 'aconB' with lower bounds 0 and upper bounds according to the values determined above.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >aconVars = {</span><span style="color: rgb(170, 4, 249);">'aconA'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconB'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model = addCOBRAVariables(model,aconVars,</span><span style="color: rgb(170, 4, 249);">'lb'</span><span >,[0;0],</span><span style="color: rgb(170, 4, 249);">'ub'</span><span >,[aconAAmount;aconBAmount]);</span></span></div></div></div><div  class = 'S10'><span></span></div><h2  class = 'S2'><span></span></h2><div  class = 'S3'><span>We further need a conversion between the used amount of aconA and the potential flux through ACONTa. We also need this for ACONTb and the same for aconB.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >linkedReactions = {</span><span style="color: rgb(170, 4, 249);">'ACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'ACONTb'</span><span >};</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >enzyme = 1:numel(aconVars)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >linkedReaction = 1:numel(linkedReactions)</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >        model = addCOBRAVariables(model,{strcat(aconVars{enzyme},</span><span style="color: rgb(170, 4, 249);">'to'</span><span >,linkedReactions{linkedReaction})},</span><span style="color: rgb(170, 4, 249);">'lb'</span><span >,0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div></div><h2  class = 'S1'><span>Adding usage efficiencies</span></h2><div  class = 'S3'><span>Now, we can add the respective constraints.</span></div><div  class = 'S3'><span>To do so, we will use the function </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%addCOBRAConstraints(model, idList, d, varargin)</span></span></div></div></div><div  class = 'S10'><span>This function will add constraints to the given </span><span style=' font-weight: bold;'>model</span><span> using a list of reactions (</span><span style=' font-weight: bold;'>idList</span><span>) involved in the constraint and right hand side values for each constraint (</span><span style=' font-weight: bold;'>d</span><span>). By default, the constraint is assumed to be a limiting constraint, and the coefficients of the reactions is assumed to be one. </span></div><div  class = 'S3'><span>i.e. if you run</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >constMod = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'FBP'</span><span >,</span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >}, 5);</span></span></div></div></div><div  class = 'S10'><span>The added constraint will restrict the sum of flow through 'FBP' and 'FBA' to 5. The function has parameters which can be used either by parameter/value pairs:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >constMod = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'FBP'</span><span >,</span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >}, 5, </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[0.2,0.4],</span><span style="color: rgb(170, 4, 249);">'dsense'</span><span >,[</span><span style="color: rgb(170, 4, 249);">'L'</span><span >]);</span></span></div></div></div><div  class = 'S10'><span>or using a parameter struct:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >params = struct();</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >params.c = [0.2,0.4];</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >params.dsense = </span><span style="color: rgb(170, 4, 249);">'L'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >constMod = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'FBP'</span><span >,</span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >}, 5, params);</span></span></div></div></div><div  class = 'S10'><span>The available parameters are:</span></div><ul  class = 'S13'><li  class = 'S14'><span style=' font-weight: bold;'>'c'</span><span> - The coefficient matrix with one column per reaction id and one row per added constraint (default: 1 for each element in </span><span style=' font-weight: bold;'>idList</span><span>)</span></li><li  class = 'S14'><span style=' font-weight: bold;'>'dsense'</span><span> - the sense vector with one element per added constraint ('E' for equality, 'L' for lower than, 'G' for greater than), or one element which is used for all constraints (default: 'L').</span></li><li  class = 'S14'><span style=' font-weight: bold;'>'ConstraintID'</span><span> - a cell array of strings with one element for each added constraint (default: ConstraintXYZ, with XYZ being the position in the ctrs vector)</span></li><li  class = 'S14'><span style=' font-weight: bold;'>checkDuplicates</span><span> - Whether to check for duplicate Constraints (they don't hurt, but they don't help). Note that duplicate IDs are still not allowed.</span></li></ul><div  class = 'S3'><span>Coming back to our example, the first constraint we add is for the amount of aconA, which should be steady (i.e. not more than made available by the aconA variable). More precisely, the amount enzyme made available for ACONTa and ACONTb, should be balanced with the amount of enzyme made available by the aconA variable.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >model = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'aconAtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconAtoACONTb'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconA'</span><span >},0, </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[-1,-1,1],</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'dsense'</span><span >,</span><span style="color: rgb(170, 4, 249);">'E'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ConstraintID'</span><span >, </span><span style="color: rgb(170, 4, 249);">'aconAAmount'</span><span >);</span></span></div></div></div><div  class = 'S10'><span>The same constraint is introduced for aconB</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >model = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'aconBtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconBtoACONTb'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconB'</span><span >},0, </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[-1,-1,1],</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'dsense'</span><span >,</span><span style="color: rgb(170, 4, 249);">'E'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ConstraintID'</span><span >, </span><span style="color: rgb(170, 4, 249);">'aconBAmount'</span><span >);</span></span></div></div></div><div  class = 'S10'><span>Next, we also add the efficiencies:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >model = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'aconAtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconBtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'ACONTa'</span><span >},0, </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[aconACit,aconBCit,-1],</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'dsense'</span><span >,</span><span style="color: rgb(170, 4, 249);">'E'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ConstraintID'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ACONTaFlux'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >model = addCOBRAConstraints(model,{</span><span style="color: rgb(170, 4, 249);">'aconAtoACONTb'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconBtoACONTb'</span><span >,</span><span style="color: rgb(170, 4, 249);">'ACONTb'</span><span >},0, </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[aconAAcon,aconBAcon,-1],</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'dsense'</span><span >,</span><span style="color: rgb(170, 4, 249);">'E'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ConstraintID'</span><span >, </span><span style="color: rgb(170, 4, 249);">'ACONTbFlux'</span><span >);</span></span></div></div></div><div  class = 'S10'><span>Finally, we have a a system, in which the two aconitase reactions are competing for the available enzymes. </span></div><h2  class = 'S2'><span>Analysing the effects of the constraints</span></h2><div  class = 'S3'><span>If we compare the results of a simple FBA optimization of the original model, and the constrained model:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S5'><span style="white-space: normal"><span >orig_sol = optimizeCbModel(model_orig)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="57CDC57F" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="258" data-hashorizontaloverflow="false" data-testid="output_3" style="max-height: 269px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">orig_sol = <span class="headerElement" style="white-space: normal; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [95×1 double]
          obj: 0.8739
        rcost: [95×1 double]
         dual: [72×1 double]
        slack: [72×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0043
        basis: [1×1 struct]
            f: 0.8739
            x: [95×1 double]
            v: [95×1 double]
            w: [95×1 double]
            y: [72×1 double]
            s: [72×1 double]
</div></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >restricted_sol = optimizeCbModel(model)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="01D9A9EB" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="314" data-hashorizontaloverflow="false" data-testid="output_4" style="max-height: 325px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">restricted_sol = <span class="headerElement" style="white-space: normal; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [101×1 double]
          obj: 0.0250
        rcost: [95×1 double]
         dual: [72×1 double]
        slack: [72×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0047
        basis: [1×1 struct]
       vars_v: [6×1 double]
       vars_w: [6×1 double]
       ctrs_y: [4×1 double]
       ctrs_s: [4×1 double]
            f: 0.0250
            x: [95×1 double]
            v: [95×1 double]
            w: [95×1 double]
            y: [72×1 double]
            s: [72×1 double]
</div></div></div></div></div></div><div  class = 'S10'><span>we can easily see that the obtained objective of the modified model is lower than that of the original model.</span></div><div  class = 'S3'><span>Values for additional variables are stored in the solution outputs </span><span style=' font-weight: bold;'>vars_v </span><span>(value used in the solution) and </span><span style=' font-weight: bold;'>vars_w </span><span>(reduced cost of the variable). For constraints the respective fields are: </span><span style=' font-weight: bold;'>ctrs_y </span><span>( fordual values for constraints) and </span><span style=' font-weight: bold;'>ctrs_s</span><span> (slacks for the constraints)</span></div><h2  class = 'S2'><span>Modifying variables and constraints</span></h2><div  class = 'S3'><span>Variables and constraints can be altered by the functions changeCOBRAVariable and changeCOBRAConstraint, respectively. Modifications on a variable include the adjustment of upper and lower bounds (by the parameters lb/ub) as well as the objective value (c) and the descriptive name (parameter 'Name').</span></div><div  class = 'S3'><span>Modifying constraints allows the adjustment of the coefficients of the constraint (parameter 'c'), the directionality (parameter 'dsense') and right hand side of the constraint ( parameter 'd') as well as the name (same as for variables). There are several ways to update the coefficients: Either a full row (including both values for the C and the E matrix - if any) needs to be provided, or the IDs along with the coefficients (similar to the situation when generating the constraints) need to be provided. It is important to note that ALL coefficients will be reset if a constraint is modified using the changeCOBRAConstraint function.</span></div><div  class = 'S3'><span>An example of these functions is provided below:</span></div><div  class = 'S3'><span>If we increase the abundance of aconA (by increasing its upper bound)</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >model = changeCOBRAVariable(model,</span><span style="color: rgb(170, 4, 249);">'aconA'</span><span >,</span><span style="color: rgb(170, 4, 249);">'ub'</span><span >,aconAAmount*2);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S16'><span style="white-space: normal"><span >less_restricted_sol = optimizeCbModel(model)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="623C8CD6" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="244" data-hashorizontaloverflow="false" data-testid="output_5" style="max-height: 261px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">less_restricted_sol = <span class="headerElement" style="white-space: normal; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [101×1 double]
          obj: 0.0254
        rcost: [101×1 double]
         dual: [76×1 double]
        slack: [76×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0030
        basis: [1×1 struct]
            x: [95×1 double]
            f: 0.0254
            y: [72×1 double]
            w: [95×1 double]
            v: [95×1 double]
</div></div></div></div></div></div><div  class = 'S10'><span>We can see, that the objective increases as more flux through the TCA is possible.</span></div><div  class = 'S3'><span>Similarily, if we reduce the efficiency of aconB on Citrate by 50%:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >model = changeCOBRAConstraints(model,</span><span style="color: rgb(170, 4, 249);">'ACONTaFlux'</span><span >,</span><span style="color: rgb(170, 4, 249);">'idList'</span><span >,{</span><span style="color: rgb(170, 4, 249);">'aconAtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'aconBtoACONTa'</span><span >,</span><span style="color: rgb(170, 4, 249);">'ACONTa'</span><span >},</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >    </span><span style="color: rgb(170, 4, 249);">'c'</span><span >,[aconACit,aconBCit*0.5,-1]);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S16'><span style="white-space: normal"><span >less_efficient_aconA = optimizeCbModel(model)</span></span></div><div  class = 'S6'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="FB7ABAB7" data-scroll-top="null" data-scroll-left="null" data-width="428" data-height="244" data-hashorizontaloverflow="false" data-testid="output_6" style="max-height: 261px; width: 458px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">less_efficient_aconA = <span class="headerElement" style="white-space: normal; font-style: italic; color: rgb(179, 179, 179); font-size: 12px;">struct with fields:</span></span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">         full: [101×1 double]
          obj: 0.0159
        rcost: [101×1 double]
         dual: [76×1 double]
        slack: [76×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0030
        basis: [1×1 struct]
            x: [95×1 double]
            f: 0.0159
            y: [72×1 double]
            w: [95×1 double]
            v: [95×1 double]
</div></div></div></div></div></div><div  class = 'S10'><span>We can again see the objective drop.</span></div><div  class = 'S3'><span>While we only used a very simple example for this tutorial, this interplay can improve predictive qualities substantially (for more have a look at e.g. [4])</span></div><div  class = 'S3'><span></span></div><h2  class = 'S2'><span>References</span></h2><div  class = 'S3'><span>[1] Thiele I, Fleming RM, Bordbar A, Schellenberger J, Palsson BØ. Functional characterization of alternate optimal solutions of Escherichia coli's transcriptional and translational machinery. Biophys J. 98(10):2072-81 (2010). </span></div><div  class = 'S3'><span>[2] The UniProt Consortium, UniProt: the universal protein knowledgebase, Nucleic Acids Res. 45: D158-D169 (2017)</span></div><div  class = 'S3'><span>[3]Jacek R. Wiśniewski, Dariusz Rakus, Quantitative analysis of the Escherichia coli proteome, Data in Brief 1, 7-11, (2014)</span></div><div  class = 'S3'><span>[4] Sánchez et al, Improving the phenotype predictions of a yeast genome-scale metabolic model by incorporating enzymatic constraints, Mol Sys Biol, 13:935 (2017)</span></div><div  class = 'S3'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Adding more complex constraints on a model 
%% Authors: Thomas Pfau,Life Sciences Reasearch Unit, University of Luxembourg, Luxembourg.
%% Reviewers: 
%% INTRODUCTION
% The COBRA Toolbox offers the possibility to add additional constraints to 
% a model that are not direct flux constraints. One example is already indicated 
% in the constraining models tutorial which introduces coupling constraints[1], 
% but there are plenty of other possibilities that can be achieved using constraints.
% 
% The tools introduced in this tutorial are mainly aimed at developers who want 
% to implement complex algorithms or formalism within the COBRA Toolbox environment, 
% but the basics are also useful to integrate specific properties of a model or 
% specific literature data. The section "The COBRA Model structure" is aimed at 
% developers and can be skipped by users who don't to go into the details.
% 
% 
%% The COBRA Model structure
% In general, the COBRA toolbox represents a model in a struct using multiple 
% fields to represent the different properties of the model. Of particular interest 
% for analysis is the stoichiometric matrix S (which represents the metabolic 
% stoichiometries of the reactions in the model). Under normal circumstances, 
% this matrix can directly be translated into the linear problem matrix A of a 
% problem struct as used by the Toolbox. However, since there are several algorithms 
% which aim at manipulating reactions or metabolites, it is important to draw 
% a distinction between the stoichiometric matrix S and any other constraints 
% or variables which can be added to the model.
% 
% Therefore 3 additional matrices exist, where such information can be stored:
% 
% E - The matrix indicating the influence of additional variables on metabolite 
% levels.
% 
% C - The matrix of additional constraints on the reactions, which are not derived 
% from the steady state condition and
% 
% D - The matrix of interactions between additional Constraints and additional 
% variables.
% 
% Overall A COBRA model will be translated into a Linear problem by combining 
% these matrices in the following way:
% 
% $$\textrm{LPproblem}\ldotp A=\;\left\lbrack \begin{array}{cc}\textrm{model}\ldotp 
% S & \textrm{model}\ldotp E\\\textrm{model}\ldotp C & \textrm{model}\ldotp D\end{array}\right\rbrack$$ 
% 
% The Cobra toolbox allows one sided inequality constraints and equality constraints 
% on a model. i.e. there is currently no mechanism to add a constraint like:
% 
% $$c_{low} \le v_{R1} \le c_{high}$$
% 
% directly. Instead it would be necessary to add the two constraints:
% 
% $$c_{low} \le v_{R1}$$
% 
% and
% 
% $$ v_{R1} \le c_{high}$$
% 
% which is commonly what solvers will translate these constraints into anyways. 
% 
% It is best practice to not manually alter the sizes of these matrices, but 
% to use the appropriate model manipulation functions (addReaction/addMultipleReactions, 
% addMetabolite/addMultipleMetabolites, addCOBRAVariables, addCOBRAConstraints). 
% These functions will ensure, that the model structure stays in sync, and all 
% necessary fields are updated appropriately. There is a convenience function 
% to create the corresponding LPproblem struct to the given model: buildOptProblemFromModel(). 
% This function will create and initialize all necessary fields and builds the 
% corresponding LP problem. Any additional modifications of the LP problem imposed 
% by an algorithm should be done on this LP and not on the original model structure.
% 
% 
% 
% 
%% PROCEDURE
% Initially, we will load a model that we want to modify by adding a few additional 
% constraints. The model used will be the simple E.Coli core model:

initCobraToolbox(false)
model = getDistributedModel('ecoli_core_model.mat');
%Create a copy for comparisons
model_orig = model;
%% 
% We will add a restriction on the activity of the two aconitase proteins present 
% in E.Coli (aconA and aconB), which catalyse two steps in the citric acid cycle 
% (see below)
% 
% 
% 
% From UNIPROT[2], the Km value of aconitase acting on citrate is 6.13 and 23.8 
% umol/min/mg for aconA and aconB respectively. For cis-aconitate the values are 
% 14.5 and 39.1 umol/min/mg, respectively. Since most models assume fluxes having 
% the unit mmol/gDW/hr we will have to convert the units to mmol/hr instead of 
% umol/min:

aconACit = 6.13 / 1000 * 60;
aconAAcon = 14.5 / 1000 * 60;
aconBCit = 23.8 / 1000 * 60;
aconBAcon = 39.1 / 1000 * 60;
%% 
% From Wiśniewski and Rakus[3] the amount of aconA and aconB per mg E.Coli sample 
% is ~4.05 pmol/mg and 95.95 pmol/mg respectively, with a weight of 97.676 kDa 
% and 93.497 kDa, respectively.

aconAmol_per_g = 4.05 * 1000 * 1e-12;
aconBmol_per_g = 95.95 * 1000 * 1e-12;
aconA_molWeight = 97.676 * 1e3;
aconB_molWeight = 93.497 * 1e3;
aconAAmount = aconAmol_per_g * aconA_molWeight / 0.3; % divided by 0.3 to account for 
                                                      % the non water fraction assuming
                                                      % 70% water.
aconBAmount = aconBmol_per_g * aconB_molWeight / 0.3;
%% 
% Now, there are two genes which code for aconitase in the ecoli core model: 
% b0118 (aconB) and b1276 (aconA). Both the citrate hydratase and the aconitatde 
% dehydratase have the same GPR rule: (b0118 or b1276), so they can both use either 
% enzyme. 

aconAgene = 'b1276';
aconBgene = 'b0118';
%% 
% We would like to add a constraint that not only restricts the activity of 
% these two reactions, but also ensures that the turnover rates are considered. 
% One mg of aconA can support a total flux of ~0.367 mmol/h of ACONTa

printRxnFormula(model,'rxnAbbrList',{'ACONTa'},'gprFlag', true);
%% 
% OR .087 mmol/hr through ACONTb.

printRxnFormula(model,'rxnAbbrList',{'ACONTb'},'gprFlag', true);
%% 
% However, it will not be able to do both at the same time. Therefore, we can 
% assume that in addition to its normal metabolites, the reaction also consumes 
% some of the enzyme (for this time step). However, both reactions can also be 
% catalysed by aconB, so they might actually not use any of aconA. Therefore, 
% we need additional constraints that represent the activity through these two 
% reactions. In addition, we need variables that represent the efficiency of one 
% mg of protein to catalyse the reactions.
% 
% Essentially, we have to encode:
% 
% A: The availability of the given proteins 
% 
% B: The efficiency of a protein to catalyse a specific reaction
% 
% C: The usage of that protein by the respective reaction.
% 
% And we will do so in this tutorial
%% Adding availability variables.
% The rxns field is intended to only represent reactions from the model, so 
% if an additional variable is required for a specific task, this variable should 
% not be generated in the rxns field, but in a distinct field for this kind of 
% variables. Therefore, we will add availability variables, which can be thought 
% of as exchange reaction for the enzymes. To do so, we use the following function

%addCOBRAVariables(model, idList, varargin)
%% 
% The function takes an existing model, and a set of variable ids (*idList*) 
% and generates all required fields for these variables. You can modify the properties 
% of the variables by the following parameters which can either be supplied as 
% parameter/value pairs:

addCOBRAVariables(model,{'newVar'},'lb', -5,'ub',3);
%% 
% or as a parameter struct:

params = struct();
params.lb = -5;
params.ub = 3;
addCOBRAVariables(model,{'newVar'},params);
%% 
% The available parameters are:
%% 
% * *'lb'* - The lower bound(s) of the variable(s) (default: -1000)
% * *'ub'* - The upper bound(s) of the variable(s) (default: 1000)
% * *'c'* - The objective coefficient(s) of the variable(s) (default: 0)
% * *'Names'* - Descriptive name(s) of the variable(s) (default: idList)
%% 
% We will now add the variables 'aconA' and 'aconB' with lower bounds 0 and 
% upper bounds according to the values determined above.

aconVars = {'aconA','aconB'};
model = addCOBRAVariables(model,aconVars,'lb',[0;0],'ub',[aconAAmount;aconBAmount]);
%% 
% 
%% 
% We further need a conversion between the used amount of aconA and the potential 
% flux through ACONTa. We also need this for ACONTb and the same for aconB.

linkedReactions = {'ACONTa','ACONTb'};
for enzyme = 1:numel(aconVars)
    for linkedReaction = 1:numel(linkedReactions)
        model = addCOBRAVariables(model,{strcat(aconVars{enzyme},'to',linkedReactions{linkedReaction})},'lb',0);
    end
end

%% Adding usage efficiencies
% Now, we can add the respective constraints.
% 
% To do so, we will use the function 

%addCOBRAConstraints(model, idList, d, varargin)
%% 
% This function will add constraints to the given *model* using a list of reactions 
% (*idList*) involved in the constraint and right hand side values for each constraint 
% (*d*). By default, the constraint is assumed to be a limiting constraint, and 
% the coefficients of the reactions is assumed to be one. 
% 
% i.e. if you run

constMod = addCOBRAConstraints(model,{'FBP','FBA'}, 5);
%% 
% The added constraint will restrict the sum of flow through 'FBP' and 'FBA' 
% to 5. The function has parameters which can be used either by parameter/value 
% pairs:

constMod = addCOBRAConstraints(model,{'FBP','FBA'}, 5, 'c',[0.2,0.4],'dsense',['L']);
%% 
% or using a parameter struct:

params = struct();
params.c = [0.2,0.4];
params.dsense = 'L';
constMod = addCOBRAConstraints(model,{'FBP','FBA'}, 5, params);
%% 
% The available parameters are:
%% 
% * *'c'* - The coefficient matrix with one column per reaction id and one row 
% per added constraint (default: 1 for each element in *idList*)
% * *'dsense'* - the sense vector with one element per added constraint ('E' 
% for equality, 'L' for lower than, 'G' for greater than), or one element which 
% is used for all constraints (default: 'L').
% * *'ConstraintID'* - a cell array of strings with one element for each added 
% constraint (default: ConstraintXYZ, with XYZ being the position in the ctrs 
% vector)
% * *checkDuplicates* - Whether to check for duplicate Constraints (they don't 
% hurt, but they don't help). Note that duplicate IDs are still not allowed.
%% 
% Coming back to our example, the first constraint we add is for the amount 
% of aconA, which should be steady (i.e. not more than made available by the aconA 
% variable). More precisely, the amount enzyme made available for ACONTa and ACONTb, 
% should be balanced with the amount of enzyme made available by the aconA variable.

model = addCOBRAConstraints(model,{'aconAtoACONTa','aconAtoACONTb','aconA'},0, 'c',[-1,-1,1],...
    'dsense','E', 'ConstraintID', 'aconAAmount');
%% 
% The same constraint is introduced for aconB

model = addCOBRAConstraints(model,{'aconBtoACONTa','aconBtoACONTb','aconB'},0, 'c',[-1,-1,1],...
    'dsense','E', 'ConstraintID', 'aconBAmount');
%% 
% Next, we also add the efficiencies:

model = addCOBRAConstraints(model,{'aconAtoACONTa','aconBtoACONTa','ACONTa'},0, 'c',[aconACit,aconBCit,-1],...
    'dsense','E', 'ConstraintID', 'ACONTaFlux');
model = addCOBRAConstraints(model,{'aconAtoACONTb','aconBtoACONTb','ACONTb'},0, 'c',[aconAAcon,aconBAcon,-1],...
    'dsense','E', 'ConstraintID', 'ACONTbFlux');
%% 
% Finally, we have a a system, in which the two aconitase reactions are competing 
% for the available enzymes. 
%% Analysing the effects of the constraints
% If we compare the results of a simple FBA optimization of the original model, 
% and the constrained model:

orig_sol = optimizeCbModel(model_orig)
restricted_sol = optimizeCbModel(model)
%% 
% we can easily see that the obtained objective of the modified model is lower 
% than that of the original model.
% 
% Values for additional variables are stored in the solution outputs *vars_v* 
% (value used in the solution) and *vars_w* (reduced cost of the variable). For 
% constraints the respective fields are: *ctrs_y* ( fordual values for constraints) 
% and *ctrs_s* (slacks for the constraints)
%% Modifying variables and constraints
% Variables and constraints can be altered by the functions changeCOBRAVariable 
% and changeCOBRAConstraint, respectively. Modifications on a variable include 
% the adjustment of upper and lower bounds (by the parameters lb/ub) as well as 
% the objective value (c) and the descriptive name (parameter 'Name').
% 
% Modifying constraints allows the adjustment of the coefficients of the constraint 
% (parameter 'c'), the directionality (parameter 'dsense') and right hand side 
% of the constraint ( parameter 'd') as well as the name (same as for variables). 
% There are several ways to update the coefficients: Either a full row (including 
% both values for the C and the E matrix - if any) needs to be provided, or the 
% IDs along with the coefficients (similar to the situation when generating the 
% constraints) need to be provided. It is important to note that ALL coefficients 
% will be reset if a constraint is modified using the changeCOBRAConstraint function.
% 
% An example of these functions is provided below:
% 
% If we increase the abundance of aconA (by increasing its upper bound)

model = changeCOBRAVariable(model,'aconA','ub',aconAAmount*2);
less_restricted_sol = optimizeCbModel(model)
%% 
% We can see, that the objective increases as more flux through the TCA is possible.
% 
% Similarily, if we reduce the efficiency of aconB on Citrate by 50%:

model = changeCOBRAConstraints(model,'ACONTaFlux','idList',{'aconAtoACONTa','aconBtoACONTa','ACONTa'},...
    'c',[aconACit,aconBCit*0.5,-1]);
less_efficient_aconA = optimizeCbModel(model)
%% 
% We can again see the objective drop.
% 
% While we only used a very simple example for this tutorial, this interplay 
% can improve predictive qualities substantially (for more have a look at e.g. 
% [4])
% 
% 
%% References
% [1] Thiele I, Fleming RM, Bordbar A, Schellenberger J, Palsson BØ. Functional 
% characterization of alternate optimal solutions of Escherichia coli's transcriptional 
% and translational machinery. Biophys J. 98(10):2072-81 (2010). 
% 
% [2] The UniProt Consortium, UniProt: the universal protein knowledgebase, 
% Nucleic Acids Res. 45: D158-D169 (2017)
% 
% [3]Jacek R. Wiśniewski, Dariusz Rakus, Quantitative analysis of the Escherichia 
% coli proteome, Data in Brief 1, 7-11, (2014)
% 
% [4] Sánchez et al, Improving the phenotype predictions of a yeast genome-scale 
% metabolic model by incorporating enzymatic constraints, Mol Sys Biol, 13:935 
% (2017)
% 
%
##### SOURCE END #####
-->
</div></body></html>